try this

what’s going on?

# For general stuff:
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.5.2
## ── Attaching packages ───────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## Warning: package 'ggplot2' was built under R version 3.5.2
## Warning: package 'tibble' was built under R version 3.5.2
## Warning: package 'tidyr' was built under R version 3.5.2
## Warning: package 'purrr' was built under R version 3.5.2
## Warning: package 'dplyr' was built under R version 3.5.2
## Warning: package 'stringr' was built under R version 3.5.2
## Warning: package 'forcats' was built under R version 3.5.2
## ── Conflicts ──────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(janitor)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library(here)
## here() starts at /Users/juansilva/Google Drive/WINTER 2020 BREN/ESM244 ADV DATA/my-site-test
## 
## Attaching package: 'here'
## The following object is masked from 'package:lubridate':
## 
##     here
library(paletteer)

# For ts stuff: 
library(tsibble)
## Warning: package 'tsibble' was built under R version 3.5.2
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:lubridate':
## 
##     interval, new_interval
## The following object is masked from 'package:dplyr':
## 
##     id
library(fable)
## Warning: package 'fable' was built under R version 3.5.2
## Loading required package: fabletools
## Warning: package 'fabletools' was built under R version 3.5.2
library(fabletools)
library(feasts)
## Warning: package 'feasts' was built under R version 3.5.2
library(forecast)
## Warning: package 'forecast' was built under R version 3.5.2
## 
## Attaching package: 'forecast'
## The following objects are masked from 'package:fabletools':
## 
##     GeomForecast, StatForecast
# For spatial stuff: 
library(sf)
## Warning: package 'sf' was built under R version 3.5.2
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
library(tmap)
## Warning: package 'tmap' was built under R version 3.5.2
library(mapview)

Monthly US energy consumption (renewables)

us_renew <- read_csv(here("data", "renewables_cons_prod.csv")) %>% 
  clean_names()
## Parsed with column specification:
## cols(
##   MSN = col_character(),
##   YYYYMM = col_double(),
##   Value = col_character(),
##   Column_Order = col_double(),
##   Description = col_character(),
##   Unit = col_character()
## )
unique(us_renew$description)
##  [1] "Wood Energy Production"             "Biofuels Production"               
##  [3] "Total Biomass Energy Production"    "Total Renewable Energy Production" 
##  [5] "Hydroelectric Power Consumption"    "Geothermal Energy Consumption"     
##  [7] "Solar Energy Consumption"           "Wind Energy Consumption"           
##  [9] "Wood Energy Consumption"            "Waste Energy Consumption"          
## [11] "Biofuels Consumption"               "Total Biomass Energy Consumption"  
## [13] "Total Renewable Energy Consumption"

We’ll focus on consumption data.

Clean up data

  • Convert description to all lowercase
  • Only keep observations for “consumption” in “description” variable
  • Remove any “total” observations from “description” variable
renew_clean <- us_renew %>% 
  mutate(description = str_to_lower(description)) %>% 
  filter(str_detect(description, pattern = "consumption")) %>% 
  filter(!str_detect(description, pattern = "total"))

Convert ‘yyyymm’ column to a date

renew_date <- renew_clean %>% 
  mutate(yr_mo_day = lubridate::parse_date_time(yyyymm, "ym")) %>% 
  mutate(month_sep = yearmonth(yr_mo_day)) %>% #coerce to tsibble `yearmonth` format
  mutate(value = as.numeric(value)) %>% 
  drop_na(month_sep, value)
## Warning: 490 failed to parse.
## Warning: NAs introduced by coercion
# Want to parse the year and month? We may use this later...
renew_parsed <-renew_date %>% 
  mutate(month = month(yr_mo_day, label = TRUE)) %>% 
  mutate(year = year(yr_mo_day))

Look at it:

renew_gg <- ggplot(data = renew_date, aes(x = month_sep, 
                                          y = value,
                                          group = description)) +
  geom_line(aes(color = description))

renew_gg

updating my colors with paleteer

renew_gg +
  scale_color_paletteer_d("ghibli::LaputaMedium")

Coerce renew_parsed to a tsibble

renew_ts <- as_tsibble(renew_parsed, key = description, index = month_sep)

Lets look at time series data in a couple diferent ways

renew_ts %>% autoplot(value)

renew_ts %>% gg_subseries(value)

# renew_ts %>% gg_season(value)

ggplot(data = renew_parsed, aes(x = month, y = value, group = year)) + 
  geom_line(aes(color = year)) +
  facet_wrap(~description,
             ncol = 1,
             scales = "free",
             strip.position = "right")

let’s look at hydroelectric consumption

hydro_ts <- renew_ts %>% 
  filter(description == "hydroelectric power consumption")

# Explore: 
hydro_ts %>% autoplot(value)

hydro_ts %>% gg_subseries(value)

# hydro_ts %>% gg_season(value)

# OK, what if gg_season() doesn't work?
# It's just a function that uses ggplot() to do things we already know how to do in ggplot()!

ggplot(hydro_ts, aes(x = month, y = value, group = year)) +
  geom_line(aes(color = year))

Calculate summary data by time using index_by()

What if we want to calculate consumption by quarter? We’ll use index_by() to tell R which “windows” to calculate a value with in.

Quarterly:

hydro_quarterly <- hydro_ts %>% 
  index_by(year_qu = ~ yearquarter(.)) %>% # monthly aggregates
  summarise(
    avg_consumption = mean(value)
  )

head(hydro_quarterly)
## # A tsibble: 6 x 2 [1Q]
##   year_qu avg_consumption
##     <qtr>           <dbl>
## 1 1973 Q1            261.
## 2 1973 Q2            255.
## 3 1973 Q3            212.
## 4 1973 Q4            225.
## 5 1974 Q1            292.
## 6 1974 Q2            290.

Decompose the hydro consumption ts data

First, let’s check the decomposition (STL):

# Find STL decomposition
dcmp <- hydro_ts %>%
  model(STL(value ~ season(window = 5)))

# View the components
# components(dcmp)

# Visualize the decomposed components
components(dcmp) %>% autoplot() +
  theme_minimal()

# Let's check out the residuals:
hist(components(dcmp)$remainder)

Explore the ACF

hydro_ts %>% 
  ACF(value) %>% 
  autoplot()

DANGER DANGWR DO A LOT OF READING

MODEL FOR FORECASTING

hydro_model <- hydro_ts %>%
  model(
    arima = ARIMA(value),
    ets = ETS(value)
  ) %>%
  fabletools::forecast(h = "2 years")

hydro_model %>% 
  autoplot(filter(hydro_ts, 
                  year(month_sep) > 2010), 
           level = NULL)

lets make a world map

# Get spatial data: 
world <- read_sf(dsn = here("data","TM_WORLD_BORDERS_SIMPL-0.3-1"), layer = "TM_WORLD_BORDERS_SIMPL-0.3") %>% clean_names()

# Quick & easy option to see those polygons (also for points, lines!)
mapview(world)
# ggplot (static)
world_base <- ggplot(data = world) +
  geom_sf(aes(fill = pop2005),
          color = NA) + 
  scale_fill_paletteer_c("viridis::viridis") +
  theme_minimal()

world_base

# Let's crop it: 
world_base +
  coord_sf(xlim = c(-20, 50), ylim = c(-40, 40), expand = FALSE)